#ifndef __eval_thresh_H__
#define __eval_thresh_H__

#include "Image.H"

// Return minimum and maximum of the input image
template <class T>
void getMinMaxImage(const pip::Image<T> &in, T &min, T &max, const pip::Image<unsigned char> &mask = pip::Image<unsigned char>())
{
  if (in.empty()) {   
    std::cerr << "No input to getMinMaxImage function!" << endl;
    return;
  }

  min = std::numeric_limits<T>::max();
  max = std::numeric_limits<T>::min();

  for (int z = 0; z < in.dim(2); ++z) {
    for (int y = 0; y < in.dim(1); ++y) {
      for (int x = 0; x < in.dim(0); ++x) {
        if(mask.empty() || mask(x, y, z)) {
          T value(in(x, y, z));

          if (value > max) {
            max = value;
          }

          if (value < min) {
            min = value;
          }
        }
      }
    }
  }
}


// evaluate the input to reference
template <typename T, typename U>
int evaluate(const pip::Image<T> &reference, 
             const pip::Image<U> &result, 
             const float &threshold = 0.5)
{
  // Check whether size of images matches
  if (reference.dim(0) != result.dim(0) || reference.dim(1) != result.dim(1) || reference.dim(2) != result.dim(2)) {   
    std::cout << "Images are of unequal size!!!" << endl;
    return 1;
  }

  // Declare true positive, true negative, false positive and false negative
  float TP(0), TN(0), FP(0), FN(0), accuracy(0),sensitivity(0), specificity(0), overlap(0); 

  // Loop over image and count TP, TN, FP, FN
  for (int z(0); z < reference.dim(2); ++z) {
    for (int y(0); y < reference.dim(1); ++y) {
      for (int x(0); x < reference.dim(0); ++x) { 
        int res((result(x, y, z) > threshold) ? 1 : 0),
          ref((int)reference(x, y, z));

        if (res == 1 && ref == 1) {
          TP++;
        } else if (res == 1 && ref == 0) {
          FP++;
        } else if (res == 0 && ref == 1) {
          FN++;
        } else if (res == 0 && ref == 0) {
          TN++;
        }
      }
    }
  }

  // Calculate outcome measures
  // Check for 0 division
  accuracy = (TP + TN) / (!(TP + TN + FP + FN) ? 1 : (TP + TN + FP + FN));
  sensitivity = TP / (!(TP + FN) ? 1 : (TP + FN));
  specificity = TN / (!(TN + FP) ? 1 : (TN + FP));
  overlap = TP / (!(TP + FP + FN) ? 1 : (TP + FP + FN));
  std::cout << "|  " << accuracy << "|  " << sensitivity<< "|  " <<specificity << "|  " << overlap << "|  " ;


  float diff(0);
  for (int y = 0; y < reference.dim(1); ++y) {
    for (int x = 0; x < reference.dim(0); ++x) {
      diff += (reference(x,y) - result(x,y)) * (reference(x,y) - result(x,y));
    }
  }

  diff = diff / (result.dim(0) * result.dim(1));

  std::cout << diff  << "sd|\n"  ;

  return 0;
}

void thresholding(const pip::Image<short> &input,
                  pip::Image<unsigned char> &output,
                  const float &threshold,
                  const pip::Image<unsigned char> &mask = pip::Image<unsigned char>());

/**
 * Performs optimal thresholding. Computes a segmentation of the input image based
 * on an optimal thresholding approach, as explained in Sonka, algorithm 6.2. A
 * separate mask can be provided to indicate the region to be segmented.
 *
 * @param input Input image
 * @param output Output image
 * @param threshold Threshold value
 * @param mask Mask image
 */
void optimal_thresholding(const pip::Image<short> &input,
                          pip::Image<unsigned char> &output,
                          float &threshold,
                          const pip::Image<unsigned char> &mask = pip::Image<unsigned char>());


float averageSquaredDifference(const pip::Image<float> &input, const pip::Image<unsigned char> &secondInput)
{
  float diff(0);
  for (int y = 0; y < input.dim(1); ++y) {
    for (int x = 0; x < input.dim(0); ++x) {
      diff += (unsigned char(input(x,y)) - secondInput(x,y)) * (unsigned char(input(x,y)) - secondInput(x,y));
    }
  }

  return (diff / (input.dim(0) * input.dim(1)));
}


#endif